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ABSTRACT 



This study concerns the fast and accurate solution of the line radiation transfer problem, under non-LTE conditions. We propose 
and evaluate an alternative iterative scheme to the classical ALI-Jacobi method, and to the more recently proposed Gauss-Seidel and 
Successive Over-Relaxation (GS/SOR) schemes. Our study is indeed based on the application of a preconditioned bi-conjugate gradi- 
ent method (BiCG-P). Standard tests, in ID plane parallel geometry and in the frame of the two-level atom model, with monochromatic 
scattering, are discussed. Rates of convergence between the previously mentioned iterative schemes are compared, as well as their 
respective timing properties. The smoothing capability of the BiCG-P method is also demonstrated. 
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1. Introduction 

The solution of the radiative transfer equation, under non-LTE 
conditions, is a classical problem in astrophysics. Many numeri- 
cal methods have been used since the beginning of the computer 
era, from difference equations methods (e.g., Feautrier 1964, 
Cuny 1967, Auer & Mihalas 1969) to iterative schemes (e.g., 
Cannon 1973, Scharmer 1981, Olson et al. 1986). 

Since the seminal paper of Olson et al. (1986), the so-called 
ALI or Accelerated K-Iteration scheme is nowadays one of 
the most popular method for solving complex radiation transfer 
problems. Using as approximate operator the diagonal of the A 
operator (see e.g., Mihalas 1978), it is a Jacobi iterative scheme. 
It has been generalized for multilevel atom problems (Rybicki & 
Hummer 1991), multi-dimensional geometries (Auer & Paletou 
1994, Auer et al. 1994, van Noort et al. 2002), and polarized 
radiation transfer (e.g., Faurobert et al. 1997, Trujillo Bueno & 
Manso Sainz 1999). 

Despite the many success of the ALI method, Trujillo Bueno 
& Fabiani Bendicho (1995) proposed a novel iterative scheme 
for the solution of non-LTE radiation transfer. They adapted 
Gauss-Seidel and Successive Over-Relaxation (GS/SOR) itera- 
tive schemes, well-known in applied mathematics as being su- 
perior, in terms of convergence rate, to the ALI-Jacobi iterative 
scheme. It is interesting to note that, besides their application 
to radiation transfer, both Jacobi and GS/SOR iterative schemes 
were proposed during the xix' /l century. 

Conjugate gradient-like, hereafter CG-like, iterative methods 
were proposed more recently by Hestenes & Stiefel (1952). They 
are of a distinct nature than ALI-Jacobi and GS/SOR methods. 
Unlike the later, CG-like methods are so-called non-stationary 
iterative methods (see e.g., Saad 2008). Very few articles dis- 
cussed its application to the radiation transfer equation (see e.g., 
Amosov & Dmitriev 2005) and, to the best of our knowledge, 
they have not been properly evaluated, so far, for astrophysical 
purposes. 



Such an evaluation is the scope of the present research note. 
This alternative approach is relevant to the quest for ever faster 
and accurate numerical methods for the solution of the non-LTE 
line radiation transfer problem. 

2. The iterative scheme 

In the two-level atom case, and with complete frequency redis- 
tribution, the non-LTE line source function is usually written as 



S(t) = (1-s)J(t) + sB(t), 



(1) 



where r is the optical depth, s is the collisional destruction prob- 
abilityfl B is the Planck function and / is the usual mean inten- 
sity: 
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(2) 



where the optical depth dependence of the specific intensity 7 v n 
(and, possibly, of the line absorption profile <p v ) is omitted for 
simplicity. 

The mean intensity is usually written as the formal solution 
of the radiative transfer equation 



/ = MS], 



(3) 



or, in other terms, the solution of the radiation transfer equation 
for a known source function (see e.g., Mihalas 1978). 

Let us define A = Id - ( 1 - s) A, the unknown x = S and the 
right-hand side b = sB. Then, the solution of the radiative trans- 
fer equation is equivalent to solving the system of equations: 



Ax = b . 



(4) 



1 The albedo a = (1 - e) is more commonly used in general studies 
of the radiation transfer equation. 
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and 




MT Zj-i = rj-i ■ 



(8) 



Fig. 1. Typical structure of the A operator, computed for the case 
of a self-emitting ID, plane parallel finite slab. Such an operator 
is, more generally, non-symmetric positive definite. For the sake 
of getting a better contrast for the off-diagonal parts, we magni- 
fied the negative values of the original operator by a factor of 5, 
and we display in fact the opposite of the modified operator. 

Among the various forms of CG-based methods, the bi- 
conjugate gradient (hereafter BiCG) method is more appropriate 
for the case of operators A which are not symmetric positive def- 
inite. It is indeed the case for our radiative transfer problem, as 
illustrated in Fig. 1 where we display the structure of the oper- 
ator A for the case of a symmetric ID grid of maximum optical 
depth r max = 2 x 10 8 , discretized with 8 spatial points per decade 
and using s = 10~ 6 . 

A general description of the BiCG gradient method can be 
found in several classical textbooks of applied mathematics (see 
e.g., Saad 2008; it is important to note that this method requires 
the use of the transpose of the A operator, A T ). However, the 
BiCG method remains efficient until the A operator becomes 
very badly conditioned i.e., for very high albedo cases, typi- 
cally when e < 10~ 6 . For the later cases, it is therefore crucial to 
switch to a more efficient scheme, using preconditioning of the 
system of equations. 

Hereafter, we derive the algorithm for the BiCG method with 
preconditioning (hereafter BiCG-P). In such a case, one seeks 
instead for a solution of: 



M~ l Ax = M~ l b, 



(5) 



where the matrix M should be easy to invert. A "natural" choice 
for the non-LTE radiation transfer problem, is to precondition 
the system of equations with the diagonal of the full operator A. 

From an initial guess for the source function, xq — S (0) , com- 
pute a residual ro such as: 



ro = b — Axi 



o ■ 



(6) 



Set the second residual r such that (ro, ro) + 0, where (a, b) 
means the inner product between vectors a and b. For all the 
cases considered by us, setting ro = ro proved to be adequate. 

Now, the BiCG-P algorithm consists in running until con- 
vergence the following iterative scheme, where j is the iterative 
scheme index. The first steps consist in solving: 



These steps are both simple and fast to compute, since M is a 
diagonal operator. Then one computes: 



Pj-i = (zJ-i'O'-i)- 



(9) 



If pj-\ = then the method fails, otherwise the algorithm contin- 
ues with the following computations. If j = 1 then, set pj = Zj-\ 
and pj = Zj-i- For j > 1, compute: 



#7-1 =Pj-l/Pj-2, 

and: 

Pj=Zj-i +Pj-iPj-i ■ 

Then make qj = Ap 7 - and qj = A T pj, as well as: 

P.H 
a; = — - . 

1 (pj>9i) 

Finally, one advances the source function according to: 
xj = Xj-i +a jPj , 
and the two residuals: 

n = n : -»,</,- 

and, 

n = o-i - a m ■ 



(10) 

(ii) 

(12) 
(13) 



(14) 



(15) 



(16) 



This process, from Eq. (0 to Eq. fl6[ , is then repeated to 
convergence. We have found that a convenient stopping criterion 
is to interrupt the iterative process when: 



^<io- 



(17) 



Mzj- 



(7) 



Indeed, this stopping criterion guarantees that the floor value of 
the true error, defined in §3.1, is always reached. 



3. Results 

The solution of the non-LTE radiative transfer problem for a 
constant property, plane -parallel ID slab is a standard test for 
any new numerical method. Indeed, the so-called monochro- 
matic scattering case allows a direct comparison of the numer- 
ical solutions to the analytical "solution", SEdd, which can be 
derived in Eddington's approximation (see e.g., the discussion 
in §5. of Chevallier et al. 2003). 
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Fig. 2. Solutions for the monochromatic scattering in an isother- 
mal atmosphere, obtained for different values of s, ranging from 
10 4 to 10~ 12 . Both the expected -\/s surface law and thermali- 
sation lengths « 1/ yfs are perfectly recovered. 



3.1. Scaling laws for non-LTE source function solutions 

In Fig. [2] we display the run of the source function with optical 
depth, for a self-emitting slab of total optical thickness r max = 2x 
10 8 , using a 9-point per decade spatial grid, a one -point angular 
quadrature with = ±1/ V3, and s values ranging from 10~ 4 to 

io- 12 . 

Both the -\/s surface law and thermalisation lengths « 1 / s/e 
are perfectly recovered, even for the numerically difficult cases 
of large (a ^ 1) albedos. 

The "true" error, that is: 



max 



s U) -s EM \ 



>Edd 



(18) 



is displayed in Fig. 3, respectively for the GS, SOR and BiCG-P 
iterative schemes. SOR was used with aw = 1.5 over-relaxation 
parameter, and s = 10~ 12 in that case. 

3.2. Timing properties 

The typical timing properties of the BiCG-P scheme, for the 
cases considered here, are summarized in Table 1. While each 
GS/SOR iteration is about 20% longer than one ALI-iteration, 
the increase of time per iteration, with the number of depth 
points, for BiCG-P is significant. However, the balance between 
computing time and reduction of the number of iterations always 
remains in favour of the BiCG-P scheme vs. GS/SOR and, a for- 
tiori, ALL 

3.3. Sensitivity to the spatial refinement 

We adopted the same slab model as the one reported in the pre- 
vious section, with s = 10~ 8 , and we allowed the grid refinement 
to vary from 5 to 20 points per decade, with a step of 3. 

Figure 4 shows the sensitivity of the method with grid re- 
finement. Increasing grid refinement corresponds to a decreasing 
value of the floor value of the true error, which demonstrates the 



Fig. 3. History of the true error vs. iteration number for e - 
10 -12 , respectively with Gauss-Seidel, SOR (with to — 1.5) and 
BiCG-P with the diagonal of A as a preconditioned 

Table 1. Timing properties of the BiCG-P iterative schemes, de- 
pending on the optical depth grid refinement expressed as the 
number of point per decade, N T , for a r ma x = 2 x 10 s , and 
e = 10~ 8 self-emitting slab. We use a time normalized to the 
time for a ALI-iteration. Note that GS/SOR is, in comparison, 
about 20% slower than ALL 
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need of a sufficient spatial resolution for the sake of accuracy on 
the numerical solutions. The number of iterative steps necessary 
to fulfill the stopping criterion defined in Eq. ( fTTb is practically 
linear in the number of depth points. 

3.4. Smoothing capability 

A strong argument in favour of the Gauss-Seidel iterative 
scheme resides in its smoothing capability. This point is partic- 
ularly important in the frame of multi-grid methods for the so- 
lution of non-LTE radiative transfer problems (see e.g., Fabiani 
Bendicho et al. 1997). 

We have reproduced the test initially proposed by Trujillo 
Bueno & Fabiani Bendicho (1995, see their Fig. 9) for the GS 
method. Here we consider a semi-infinite slab of maximum op- 
tical thickness 10 5 , sampled with a 9-point per decade grid, and 
s = 10~ 6 . The test consists in injecting ab initio a high-frequency 
component into the error, 5 (0) - 5 (oo) , on the initial guess for the 
source function. In Fig. 5, the initial error we used is easily rec- 
ognized as the hyperbolic tangent-like curve, to which a high- 
frequency component was added. 

After the very first iterate of the BiCG-P iterative process, 
the error has an amplitude of * 0.75 and, more importantly, it 
is already completely smoothed. Other error curves displayed in 
Fig. 5, of continuously decreasing amplitudes, are the successive 
errors at iterates 2 to 5. 

This simple numerical experiment demonstrates that BiCG-P 
can also be used as an efficient smoother for multi-grid methods. 
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Fig. 4. Evolution of the true error for varying spatial grid refine- 
ments, from 5 to 20 points per decade, with a step of 3. The 
value of the true error plateau decreases with an increasing grid 
refinement. 




tween such over-computing time and an improved convergence 
rate, with respect to the ones of GS/SOR iterative schemes, is 
currently being conducted. 

Multi-level and multi-dimensional cases will be considered 
further. However, the use of BiCG-P does not require the cum- 
bersome modifications of multi-dimensional formal solvers re- 
quired by GS/SOR methods (see e.g., Leger et al. 2007). 

Acknowledgements. We are grateful to Drs. Bernard Rutily, and Loic Chevallier 
for numerous discussions about the solution of the radiative transfer equation. 
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Fig. 5. Initialization and evolution of the error, S (j) - S (oo \ during 
the 5 first iterations of the BiCG-P scheme. The high-frequency 
component introduced ab initio in the source function is re- 
moved after the very first iteration of BiCG-P. Other curves, of 
continuously decreasing amplitudes, correspond to iterations 2 
to 5. 



4. Conclusion 

We propose an alternative method for the solution of the non- 
LTE radiation transfer problem. Preliminary but standards tests 
for the two-level atom case in ID plane parallel geometry are 
successful. In such a case, the timing properties of the precondi- 
tioned BiCG method are comparable to the ones of ALI, GS and 
SOR stationnary iterative methods. However the convergence 
rate of BiCG-P is signifiantly better. 

The main potential drawback of the BiCG-P method is in 
the usage of the A T operator. A detailed trade-off analysis be- 



